!*************************************************************************
subroutine task3d(vertex2D,nVer2D,nInCat2D,vertex3D,params,nVer3D,nInCat3D,&
    igT,jgT,ggT,ig,jg,ggu,ggl,di,f,l1,v1,v2,v3,v4,v5,ggl2,ggu2,di2,vnRZ,vnXYZ,&
    igT3D,jgT3D,ggT3D,l13d,va,ig3D,jg3D,igEarth,jgEarth,a,nodesforedges,a0,&
    vn2d_dat,vn3d_dat,va_dat,ax_dat,ay_dat,az_dat,a_0000,aXYZ)
    common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,m2D,mTerm2d,nodesNotTermCount2D
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,edgesNotTermCount3D,&
    edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimEarth/nodesNotTermEarth,mEarth
    common /files_common/ infite_0,num_slae,log_txt,typeTask_txt
    integer elemsCount2D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,elementsEarth
    integer typeTask,calcWK,calcVP,vn2d_dat,vn3d_dat,va_dat,log_txt,ax_dat,ay_dat,az_dat,a_0000,aXYZ,typeTask_txt
    real*8 eps2d,eps3d,eps3dNS,radiusRED,vertex2D,vertex3D,params,ggT,ggu,ggl,di,&
    v1,v2,v3,v4,v5,ggl2,ggu2,di2,vnRZ,vnXYZ,ggT3D,va,f,a,a0    
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    vertex3D(3,nodesCount3D),params(2,nMaterials),nVer3D(14,elemsCount3D),nInCat3D(elemsCount3D),&
    igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),jg(m2d),ggu(m3d),&
    ggl(m3d),di(nodesNotTermCount3D),l1(nodesBoundCount2D),v1(nodesNotTermCount3D),&
    v2(nodesNotTermCount3D),v3(nodesNotTermCount3D),v4(nodesNotTermCount3D),v5(nodesNotTermCount3D),&
    ggl2(m3d),ggu2(m3d),di2(nodesNotTermCount3D),vnRZ(nodesCount2D),vnXYZ(nodesCount3D),&
    igT3D(nodesTermCount3D+1),jgT3D(mTerm3d),ggT3D(mTerm3d),l13d(nodesBoundCount3D),va(nodesCount3D),&
    ig3D(nodesNotTermCount3D+1),jg3D(m2d),igEarth(nodesNotTermEarth+1),jgEarth(mEarth),&
    f(nodesNotTermCount3D,3),a(nodesCount3D,3),nodesforedges(3,edgesCount3D),a0(edgesCount3D)

    !calculating normal component of V in RZ and recalculating it to XYZ
    call V_n(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,di,&
    f(1,1),l1,v1,v2,v3,v4,v5,ggl2,ggu2,di2,vnRZ,vertex3D,nVer3D,nInCat3D,vnXYZ,&
    igT3D,jgT3D,ggT3D,'SP',vn2d_dat,vn3d_dat,log_txt)        

    !calculating anomal component of V
    call V_a(params,ig3D,jg3D,ggu,ggl,di,f(1,1),l13d,v1,v2,v3,v4,v5,ggl2,ggu2,di2,&
    vnXYZ,vertex3D,nVer3D,nInCat3D,igT3D,jgT3D,ggT3D,va,igEarth,jgEarth,'SP',va_dat,&
    log_txt)        

    !calculating anomal component of vector A0
    call Aa_0(params,ig3D,jg3D,ggu,ggl,di,f,l13d,v1,v2,v3,v4,v5,ggl2,ggu2,di2,vnXYZ,&
    vertex3D,nVer3D,nInCat3D,igT3D,jgT3D,ggT3D,va,a(1,1),a(1,2),a(1,3),nodesforedges,&
    a0,ax_dat,ay_dat,az_dat,a_0000,aXYZ,log_txt)

    !calculation rotAx (validation test)
!    call rotA(vertex2D,vertex3D,nVer3D,a(1,1),a(1,2),a(1,3))  
end
!*************************************************************************
subroutine V_n(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,&
    di,f,l1,v1,v2,v3,v4,v5,ggl2,ggu2,di2,vnRZ,vertex3D,nVer3D,nInCat3D,vnXYZ,&
    igT3D,jgT3D,ggT3D,prefix,vn2d_dat,vn3d_dat,log_txt)
    common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,m2D,mTerm2d,nodesNotTermCount2D
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,edgesNotTermCount3D,&
    edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    integer elemsCount2D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D
    integer elem1,elem2,typeTask,calcWK,calcVP,vn2d_dat,vn3d_dat,log_txt
    real*8 vertex2D,vertex3D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,vnRZ,vnXYZ
    real*8 eps2d,eps3d,eps3dNS,radiusRED,ggT3D
    character*2 prefix
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(2,nMaterials)
    dimension igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),f(nodesNotTermCount2D),l1(nodesBoundCount2D),&
    v1(nodesNotTermCount2D),v2(nodesNotTermCount2D),v3(nodesNotTermCount2D),v4(nodesNotTermCount2D),&
    v5(nodesNotTermCount2D),ggl2(m2d),ggu2(m2d),di2(nodesNotTermCount2D),vnRZ(nodesCount2D),&
    vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),nInCat3D(elemsCount3D),vnXYZ(nodesCount3D),&
    igT3D(nodesTermCount3D+1),jgT3D(mTerm3d),ggT3D(mTerm3d)

    call clearVectorDouble(vnXYZ,nodesCount3D)

    !Assembling global matrix
    call globalForVn(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,ggT,igT,jgT,elem1,elem2)

    !Taking into account boundary condition
    call firstCondition(di,f,l1,nodesNotTermCount2D,nodesBoundCount2D,1)
    call secondConditionForVn(f,nVer2D(1,elem1),nVer2D(3,elem1),nVer2D(2,elem2),nVer2D(4,elem2))

    !Solving SLAE
    call LOS3(vnRZ,f,ig,jg,ggl,ggu,di,nodesNotTermCount2D,m2d,&
    maxiter,eps2d,v1,v2,v3,v4,ggl2,ggu2,di2,v5,0,iter,residual)
    call calculateValueForNotTermNodes(vnRZ,igT,jgT,ggT,nodesCount2D,&
    mTerm2d,nodesTermCount2D,nodesNotTermCount2D)

    !Recalculating solution to XYZ
    call v3d(vertex3D,nVer3D,vertex2D,nVer2D,vnRZ,vnXYZ)
    call calculateValueForNotTermNodes(vnXYZ,igT3D,jgT3D,ggT3D,nodesCount3D,&
    mTerm3d,nodesTermCount3D,nodesNotTermCount3D)

    !Writing solution in file
    call writeToFileDouble(vn2d_dat,vnRZ,nodesCount2D)
    call writeToFileDouble(vn3d_dat,vnXYZ,nodesCount3D)
    close(vn2d_dat)
    close(vn3d_dat)

2   format(A,A,i5,' residual =',E)
    write(*,2) prefix,': end V_n:       iter =',iter,residual
    write(log_txt,2)prefix,': end V_n:       iter =',iter,residual
end
!*************************************************************************
subroutine V_a(params,ig,jg,ggu,ggl,di,f,l1,v1,v2,v3,v4,v5,ggl2,ggu2,di2,vnXYZ,vertex3D,&
    nVer3D,nInCat3D,igT3D,jgT3D,ggT3D,va,igEarth,jgEarth,prefix,va_dat,log_txt)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,edgesNotTermCount3D,&
    edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimEarth/nodesNotTermEarth,mEarth
    integer va_dat,log_txt
    integer edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,elementsEarth
    real*8 vertex3D,params,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,vnXYZ
    character*2 prefix
    dimension params(2,nMaterials),ig(nodesNotTermCount3D+1),jg(m3d),ggu(m3d),ggl(m3d),&
    di(nodesNotTermCount3D),f(nodesNotTermCount3D),l1(nodesBoundCount3D),v1(nodesNotTermCount3D),&
    v2(nodesNotTermCount3D),v3(nodesNotTermCount3D),v4(nodesNotTermCount3D),v5(nodesNotTermCount3D),&
    ggl2(m3d),ggu2(m3d),di2(nodesNotTermCount3D),vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),&
    nInCat3D(elemsCount3D),vnXYZ(nodesCount3D),igT3D(nodesTermCount3D+1),jgT3D(mTerm3d),ggT3D(mTerm3d),&
    va(nodesCount3D),igEarth(nodesNotTermEarth+1),jgEarth(mEarth)
    real*8 eps2d,eps3d,eps3dNS,radiusRED,ggT3D,va

    call clearVectorDouble(va,nodesCount3D)

    !Assembling global matrix
    call globalForVa(ggl,ggu,di,f,ig,jg,vertex3D,nVer3D,nInCat3D,params,ggT3D,igT3D,jgT3D,vnXYZ)

    !Taking into account boundary condition
    call firstCondition(di,f,l1,nodesNotTermCount3D,nodesBoundCount3D,1)
    k=1
    do i=1,nodesNotTermEarth
        j=ig(i)
        do while(j<=ig(i+1)-1.and.jg(j)<=nodesNotTermEarth)
            ggl2(k)=ggl(j)
            ggu2(k)=ggu(j)
            j=j+1
            k=k+1
        enddo
    enddo

    !Solving SLAE
    call LOS2(va,f,igEarth,jgEarth,ggl2,ggu2,di,nodesNotTermEarth,mEarth,&
    maxiter,eps3d,v1,v2,v3,v4,v5,0,iter,residual)

    call calculateValueForNotTermNodes(va,igT3D,jgT3D,ggT3D,nodesCount3D,&
    mTerm3d,nodesTermCount3D,nodesNotTermCount3D)

    !Writing solution in file
    call writeToFileDouble(va_dat,va,nodesCount3D)
    close(va_dat)

2   format(A,A,i5,' residual =',E)
    write(*,2) prefix,': end V_a:       iter =',iter,residual
    write(log_txt,2)prefix,': end V_a:       iter =',iter,residual
end
!*************************************************************************
subroutine Aa_0(params,ig,jg,ggu,ggl,di,f,l1,v1,v2,v3,v4,v5,ggl2,ggu2,di2,&
    vnXYZ,vertex3D,nVer3D,nInCat3D,igT3D,jgT3D,ggT3D,va,aX,aY,aZ,nodesforedges,&
    a0,ax_dat,ay_dat,az_dat,a_0000,aXYZ,log_txt)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    integer ax_dat,ay_dat,az_dat,a_0000,log_txt,aXYZ
    integer edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,calcWK,calcVP,&
    typeTask
    real*8 vertex3D,params,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,vnXYZ
    dimension params(2,nMaterials),ig(nodesNotTermCount3D+1),jg(m3d),ggu(m3d),ggl(m3d),&
    di(nodesNotTermCount3D),l1(nodesBoundCount3D),v1(nodesNotTermCount3D),&
    v2(nodesNotTermCount3D),v3(nodesNotTermCount3D),v4(nodesNotTermCount3D),v5(nodesNotTermCount3D),&
    ggl2(m3d),ggu2(m3d),di2(nodesNotTermCount3D),vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),&
    nInCat3D(elemsCount3D),vnXYZ(nodesCount3D),igT3D(nodesTermCount3D+1),jgT3D(mTerm3d),ggT3D(mTerm3d),&
    va(nodesCount3D),f(nodesNotTermCount3D,3),aX(nodesCount3D),aY(nodesCount3D),aZ(nodesCount3D),&
    nodesforedges(3,edgesCount3D),a0(edgesCount3D)
    real*8 eps2d,eps3d,eps3dNS,radiusRED,ggT3D,va,aX,aY,aZ,a0

    !Assembling global matrix
    call globalForAa_0(ggl,ggu,di,f,ig,jg,vertex3D,nVer3D,nInCat3D,params,ggT3D,igT3D,jgT3D,vnXYZ,va)

    !Taking into account boundary condition
    call firstCondition(di,f,l1,nodesNotTermCount3D,nodesBoundCount3D,3) 

    !Solving SLAE
2   format(A,i5,' residual =',E)
    call LOS2(aX,f(1,1),ig,jg,ggl,ggu,di,nodesNotTermCount3D,m3d,&
    maxiter,eps3d,v1,v2,v3,v4,v5,0,iter,residual)
    write(*,2) 'SP: end A_aX:      iter =',iter,residual
    write(log_txt,2)'SP: end A_aX:      iter =',iter,residual
    call LOS2(aY,f(1,2),ig,jg,ggl,ggu,di,nodesNotTermCount3D,m3d,&
    maxiter,eps3d,v1,v2,v3,v4,v5,0,iter,residual)
    write(*,2) 'SP: end A_aY:      iter =',iter,residual
    write(log_txt,2)'SP: end A_aY:      iter =',iter,residual
    call LOS2(aZ,f(1,3),ig,jg,ggl,ggu,di,nodesNotTermCount3D,m3d,&
    maxiter,eps3d,v1,v2,v3,v4,v5,0,iter,residual)
    write(*,2) 'SP: end A_aZ:      iter =',iter,residual
    write(log_txt,2)'SP: end A_aZ:      iter =',iter,residual

    call calculateValueForNotTermNodes(aX,igT3D,jgT3D,ggT3D,nodesCount3D,&
    mTerm3d,nodesTermCount3D,nodesNotTermCount3D)
    call calculateValueForNotTermNodes(aY,igT3D,jgT3D,ggT3D,nodesCount3D,&
    mTerm3d,nodesTermCount3D,nodesNotTermCount3D)
    call calculateValueForNotTermNodes(aZ,igT3D,jgT3D,ggT3D,nodesCount3D,&
    mTerm3d,nodesTermCount3D,nodesNotTermCount3D)
    call calculateA0(aX,aY,aZ,a0,vertex3D,nodesforedges)

    !Writing solution in file
    call writeToFileDouble(ax_dat,aX,nodesCount3D)
    call writeToFileDouble(ay_dat,aY,nodesCount3D)
    call writeToFileDouble(az_dat,aZ,nodesCount3D)
    call writeToFileDouble(a_0000,a0,edgesNotTermCount3D)
    do j=1,nodesCount3D
            write(aXYZ,rec=j)aX(j),aY(j),aZ(j)
    enddo
    close(ax_dat)
    close(ay_dat)
    close(az_dat)
    close(a_0000)
    close(aXYZ)
end
!*************************************************************************
subroutine rotA(vertex2D,vertex3D,nVer3D,aX,aY,aZ)
    common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,m2D,mTerm2d,nodesNotTermCount2D
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    integer elemsCount2D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,searchElement3D
    real*8 vertex2D,vertex3D,aX,aY,aZ,dAz,dAy,x,y,z,rotFx
    dimension vertex2D(2,nodesCount2D),vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),&
    aX(nodesCount3D),aY(nodesCount3D),aZ(nodesCount3D)

!    open(500,file='.\calculations\ax1.dat',access='direct',recl=8,mode='read')
!    call readFromFileDouble(500,aX,nodesCount3D)
!    close(500)
!    open(500,file='.\calculations\ay1.dat',access='direct',recl=8,mode='read')
!    call readFromFileDouble(500,aY,nodesCount3D)
!    close(500)
!    open(500,file='.\calculations\az1.dat',access='direct',recl=8,mode='read')
!    call readFromFileDouble(500,aZ,nodesCount3D)
!    close(500)

    open(500,file='.\calculations\rotAx.dat',access='direct',recl=8,status='replace',mode='write')
    do i=1,nodesCount2D
        x = 0.d0
        y = vertex2D(1,i)!1400.d0 !
        z = vertex2D(2,i)!-2000.d0 !
        nEl = searchElement3D(vertex3D,nVer3D,x,y,z,nodesCount3D,elemsCount3D)
        write(500,rec=i) rotFx(vertex3D,nVer3D,aY,aZ,x,y,z,nEl)
    enddo
    close(500)
end
!*************************************************************************
real*8 function rotFx(vertex3D,nVer3D,aY,aZ,x,y,z,nEl)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    integer edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,searchElement3D
    real*8 vertex3D,aY,aZ,x,y,z,rotX,ValueInPoint3D,dAz,dAy,h1,h2,k1,k2,PartitionInPoint
    dimension vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),aY(nodesCount3D),&
    aZ(nodesCount3D),rotX(8),dAz(8),dAy(8)

    do i=1,8
        rotX(i) = PartitionInPoint(aY,aZ,vertex3D,nVer3D,x,y,z,nEl)/(4.d0*3.1415926d0*1d-7)
    enddo
    rotFx = -ValueInPoint3D(rotX,vertex3D,nVer3D,x,y,z,nEl)    
end
!*************************************************************************
real*8 function PartitionInPoint(vY,vZ,vertex,nVer,x,y,z,num)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    real*8 vertex,vY,vZ,x,y,z,x1,y1,z1,x2,y2,z2,functionX,functionY,functionZ,dfunctionY,dfunctionZ,hx,hy,hz
    integer elemsCount3D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,edgesBoundCount3D
    dimension vertex(3,nodesCount3D),vY(nodesCount3D),vZ(nodesCount3D),nVer(14,elemsCount3D)

    x1 = vertex(1,nVer(1,num))
    x2 = vertex(1,nVer(2,num))
    y1 = vertex(2,nVer(1,num))
    y2 = vertex(2,nVer(3,num))
    z1 = vertex(3,nVer(1,num))
    z2 = vertex(3,nVer(5,num))

    hx=x2-x1
    hy=y2-y1
    hz=z2-z1
    PartitionInPoint=0.d0
    do i=1,8
        PartitionInPoint=PartitionInPoint+functionX(x1,x2,hx,x,i)*&
                                          dfunctionY(y1,y2,hy,y,i)*&
                                          functionZ(z1,z2,hz,z,i)*&
                                          vZ(nVer(i,num))-&
                                          functionX(x1,x2,hx,x,i)*&
                                          functionY(y1,y2,hy,y,i)*&
                                          dfunctionZ(z1,z2,hz,z,i)*&
                                          vY(nVer(i,num))
    enddo
end
!*************************************************************************
!Assembling global matrix(Vn)
subroutine globalForVn(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,ggT,&
    igT,jgT,elem1,elem2)
    common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,m2D,mTerm2d,nodesNotTermCount2D
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,edgesNotTermCount3D,&
    edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    integer elemsCount2D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,edgeTerm2d
    integer elem1,elem2,calcWK,calcVP,typeTask
    real*8 vertex2D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,vnRZ
    real*8 g,r1,r2,z1,z2,sigma,b,eps2d,eps3d,eps3dNS,radiusRED
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(2,nMaterials)
    dimension g(4),b(4,4),igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),&
    ig(nodesNotTermCount2D+1),jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),&
    f(nodesNotTermCount2D)

    call clearVectorDouble(g,4)
    call clear(ggl,ggu,di,f,nodesNotTermCount2D,m2d)
    elem1=0
    elem2=0
    do i=1,elemsCount2D
        r1=vertex2D(1,nVer2D(1,i))
        r2=vertex2D(1,nVer2D(2,i))
        z1=vertex2D(2,nVer2D(1,i))
        z2=vertex2D(2,nVer2D(3,i))
        sigma=params(2,nInCat2D(i))
        if(r1.eq.0.d0.and.z2.eq.0.d0) elem1=i
        if(r2.eq.radiusRED.and.z2.eq.0.d0) elem2=i
        if(r1.eq.0.d0)r1=1d-3

        !Assembling local matrix
        call LocalMatrixForVn(b,r1,r2,z1,z2,sigma)
        !Adding local matrix to global
        call AddLocal(di,ig,jg,ggl,ggu,f,b,g,nVer2D(1,i),nodesNotTermCount2D,&
        m2d,4,igT,jgT,ggT,nodesTermCount2D,mTerm2d,1)
    enddo
end
!*************************************************************************
!Assembling global matrix(Va)
subroutine globalForVa(ggl,ggu,di,f,ig,jg,vertex3D,nVer3D,nInCat3D,params,&
    ggT3D,igT3D,jgT3D,vnXYZ)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    integer edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,&
    elementsEarth,typeTask,calcWK,calcVP
    real*8 vertex3D,params,ggu,ggl,di,f,ggl2,ggu2,di2,vnXYZ
    real*8 ggT3D,v,b,x1,x2,y1,y2,z1,z2,g,sigma0,sigma,eps2d,eps3d,eps3dNS,radiusRED
    dimension params(2,nMaterials),ig(nodesNotTermCount3D+1),jg(m3d),ggu(m3d),ggl(m3d),&
    di(nodesNotTermCount3D),f(nodesNotTermCount3D),vertex3D(3,nodesCount3D),&
    nVer3D(14,elemsCount3D),nInCat3D(elemsCount3D),vnXYZ(nodesCount3D),&
    igT3D(nodesTermCount3D+1),jgT3D(mTerm3d),ggT3D(mTerm3d),v(8),b(8,8),g(8)

    call clear(ggl,ggu,di,f,nodesNotTermCount3D,m3D)
    do i=1,elemsCount3D
        x1=vertex3D(1,nVer3D(1,i))
        x2=vertex3D(1,nVer3D(2,i))
        y1=vertex3D(2,nVer3D(1,i))
        y2=vertex3D(2,nVer3D(3,i))
        z1=vertex3D(3,nVer3D(1,i))
        z2=vertex3D(3,nVer3D(5,i))
        sigma0=params(2,nInCat3D(i))
        sigma=params(1,nInCat3D(i))
        do j=1,8
            v(j)=vnXYZ(nVer3D(j,i))
        enddo

        !Assembling local matrix
        call LocalMatrixForVa(b,g,x1,x2,y1,y2,z1,z2,sigma0,sigma,v)

        !Adding local matrix to global
        call AddLocal(di,ig,jg,ggl,ggu,f,b,g,nVer3D(1,i),nodesNotTermCount3D,m3D,8,igT3D,jgT3D,ggT3D,&
        nodesTermCount3D,mTerm3d,1)
    enddo
end
!*************************************************************************
!Assembling global matrix(Aa_0)
subroutine globalForAa_0(ggl,ggu,di,f,ig,jg,vertex3D,nVer3D,nInCat3D,params,ggT3D,igT3D,jgT3D,vnXYZ,va)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,edgesNotTermCount3D,&
    edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    integer edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,typeTask,calcWK,calcVP
    real*8 vertex3D,params,ggu,ggl,di,f,ggl2,ggu2,di2,vnXYZ
    dimension params(2,nMaterials),ig(nodesNotTermCount3D+1),jg(m3d),ggu(m3d),ggl(m3d),&
    di(nodesNotTermCount3D),vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),nInCat3D(elemsCount3D),&
    vnXYZ(nodesCount3D),igT3D(nodesTermCount3D+1),jgT3D(mTerm3d),ggT3D(mTerm3d),v1(8),v2(8),&
    b(8,8),g(8,3),f(nodesNotTermCount3D,3),va(nodesCount3D)
    real*8 ggT3D,v1,v2,b,x1,x2,y1,y2,z1,z2,g,sigma0,sigma,f_mu,va,eps2d,eps3d,eps3dNS,radiusRED

    call clear(ggl,ggu,di,f(1,1),nodesNotTermCount3D,m3d)
    call clearVectorDouble(f(1,2),nodesNotTermCount3D)
    call clearVectorDouble(f(1,3),nodesNotTermCount3D)
    do i=1,elemsCount3D
        x1=vertex3D(1,nVer3D(1,i))
        x2=vertex3D(1,nVer3D(2,i))
        y1=vertex3D(2,nVer3D(1,i))
        y2=vertex3D(2,nVer3D(3,i))
        z1=vertex3D(3,nVer3D(1,i))
        z2=vertex3D(3,nVer3D(5,i))
        sigma0=params(2,nInCat3D(i))
        sigma=params(1,nInCat3D(i))
        f_mu=1.d0/(4.d0*3.1415926d0*1d-7)

        do j=1,8
            v1(j)=vnXYZ(nVer3D(j,i))
            v2(j)=va(nVer3D(j,i))
        enddo

        !Assembling local matrix
        call LocalMatrixVercorForAa_0(b,g(1,1),g(1,2),g(1,3),x1,x2,y1,y2,z1,z2,sigma0,sigma,f_mu,v1,v2)

        !Adding local matrix to global
        call AddLocal(di,ig,jg,ggl,ggu,f,b,g,nVer3D(1,i),nodesNotTermCount3D,m3d,8,igT3D,jgT3D,ggT3D,&
        nodesTermCount3D,mTerm3d,3)
    enddo
end
!*************************************************************************
!Assembling local matrix for Vn
subroutine LocalMatrixForVn(b,r1,r2,z1,z2,sigma)
    dimension b(4,4),br(2,2),cr(2,2),bz(2,2),cz(2,2)
    real*8 b,r1,r2,z1,z2,br,cr,bz,cz,sigma 
    call local_bc_r(br,cr,r1,r2)
    call local_bc_z(bz,cz,z1,z2)

        do i=1,4
            do j=1,4
              mi=mu(i)
              mj=mu(j)
              ni=nu(i)
              nj=nu(j)
              b(i,j)=(br(mi,mj)*cz(ni,nj)+bz(ni,nj)*cr(mi,mj))*sigma 
            enddo
        enddo
end
!*************************************************************************
!Assembling local matrix for Va
subroutine LocalMatrixForVa(b,g,x1,x2,y1,y2,z1,z2,sigma_n,sigma,v)
    dimension b1(2,2),c1(2,2),b(8,8),g(8),v(8)
    real*8 b1,c1,b,hx,hy,hz,g,x1,x2,y1,y2,z1,z2,v,sigma,sigma_n
    
    b1(1,1)=1.d0
    b1(1,2)=-1.d0
    b1(2,1)=-1.d0
    b1(2,2)=1.d0

    c1(1,1)=2.d0/6.d0
    c1(1,2)=1.d0/6.d0
    c1(2,1)=1.d0/6.d0
    c1(2,2)=2.d0/6.d0

    hx=x2-x1
    hy=y2-y1
    hz=z2-z1
    do i=1,8
        do j=1,8
            mi=mu_3d(i)
            mj=mu_3d(j)
            ni=nu_3d(i)
            nj=nu_3d(j)
            ki=ku_3d(i)
            kj=ku_3d(j)
            b(i,j)=(1.d0/hx*b1(mi,mj)*hy*c1(ni,nj)*hz*c1(ki,kj)+&
            hx*c1(mi,mj)*1.d0/hy*b1(ni,nj)*hz*c1(ki,kj)+&
            hx*c1(mi,mj)*hy*c1(ni,nj)*1.d0/hz*b1(ki,kj))
        enddo
    enddo

    call mult_sq(g,b,v,8)
    do i=1,8
        g(i)=g(i)*(sigma_n-sigma)
        do j=1,8
            b(i,j)=b(i,j)*sigma
        enddo
    enddo
end
!*************************************************************************
!Assembling local matrix for Aa_0
subroutine LocalMatrixVercorForAa_0(b,gx,gy,gz,x1,x2,y1,y2,z1,z2,sigma0,sigma,f_mu,v1,v2)
    dimension b(8,8),dx(8,8),v1(8),v2(8),gx(8),gy(8),gz(8)
    real*8 b,v1,v2,gx,gy,gz,x1,x2,y1,y2,z1,z2,sigma,sigma0,hx,hy,hz,f_mu,dx

    hx=x2-x1
    hy=y2-y1
    hz=z2-z1
    call localMatrixForAa_0(b,hx,hy,hz,f_mu)
    call local_dx(gx,dx,hy,hz,v1,v2,sigma,sigma0)
    call local_dy(gy,dx,hx,hz,v1,v2,sigma,sigma0)
    call local_dz(gz,dx,hx,hy,v1,v2,sigma,sigma0)
end
!*************************************************************************
subroutine localMatrixForAa_0(b,hx,hy,hz,f_mu)
    dimension b1(2,2),c1(2,2),b(8,8)
    real*8 b1,c1,b,hx,hy,hz,f_mu
    b1(1,1)=1.d0
    b1(1,2)=-1.d0
    b1(2,1)=-1.d0
    b1(2,2)=1.d0

    c1(1,1)=2.d0/6.d0
    c1(1,2)=1.d0/6.d0
    c1(2,1)=1.d0/6.d0
    c1(2,2)=2.d0/6.d0

    do i=1,8
        do j=1,8
            mi=mu_3d(i)
            mj=mu_3d(j)
            ni=nu_3d(i)
            nj=nu_3d(j)
            ki=ku_3d(i)
            kj=ku_3d(j)
            b(i,j)=(1.d0/hx*b1(mi,mj)*hy*c1(ni,nj)*hz*c1(ki,kj)+&
            hx*c1(mi,mj)*1.d0/hy*b1(ni,nj)*hz*c1(ki,kj)+&
            hx*c1(mi,mj)*hy*c1(ni,nj)*1.d0/hz*b1(ki,kj))*f_mu
        enddo
    enddo
end
!*************************************************************************
subroutine local_dx(g,dx,h1,h2,LocalVn,LocalVa,sigma,sigma_n)
    dimension dx(8,8),LocalVn(8),LocalVa(8),g(8),g_(8)
    real*8 dx,LocalVn,LocalVa,g,g_,sigma,sigma_n,h1,h2

    dx(1,1)=-4.d0
    dx(1,2)=4.d0
    dx(1,3)=-2.d0
    dx(1,4)=2.d0
    dx(1,5)=-2.d0
    dx(1,6)=2.d0
    dx(1,7)=-1.d0
    dx(1,8)=1.d0

    dx(2,1)=-4.d0
    dx(2,2)=4.d0
    dx(2,3)=-2.d0
    dx(2,4)=2.d0
    dx(2,5)=-2.d0
    dx(2,6)=2.d0
    dx(2,7)=-1.d0
    dx(2,8)=1.d0

    dx(3,1)=-2.d0
    dx(3,2)=2.d0
    dx(3,3)=-4.d0
    dx(3,4)=4.d0
    dx(3,5)=-1.d0
    dx(3,6)=1.d0
    dx(3,7)=-2.d0
    dx(3,8)=2.d0

    dx(4,1)=-2.d0
    dx(4,2)=2.d0
    dx(4,3)=-4.d0
    dx(4,4)=4.d0
    dx(4,5)=-1.d0
    dx(4,6)=1.d0
    dx(4,7)=-2.d0
    dx(4,8)=2.d0

    dx(5,1)=-2.d0
    dx(5,2)=2.d0
    dx(5,3)=-1.d0
    dx(5,4)=1.d0
    dx(5,5)=-4.d0
    dx(5,6)=4.d0
    dx(5,7)=-2.d0
    dx(5,8)=2.d0

    dx(6,1)=-2.d0
    dx(6,2)=2.d0
    dx(6,3)=-1.d0
    dx(6,4)=1.d0
    dx(6,5)=-4.d0
    dx(6,6)=4.d0
    dx(6,7)=-2.d0
    dx(6,8)=2.d0

    dx(7,1)=-1.d0
    dx(7,2)=1.d0
    dx(7,3)=-2.d0
    dx(7,4)=2.d0
    dx(7,5)=-2.d0
    dx(7,6)=2.d0
    dx(7,7)=-4.d0
    dx(7,8)=4.d0

    dx(8,1)=-1.d0
    dx(8,2)=1.d0
    dx(8,3)=-2.d0
    dx(8,4)=2.d0
    dx(8,5)=-2.d0
    dx(8,6)=2.d0
    dx(8,7)=-4.d0
    dx(8,8)=4.d0

    do i=1,8
        do j=1,8
            dx(i,j)=dx(i,j)/72.d0*h1*h2
        enddo
    enddo

    call mult_sq(g,dx,LocalVn,8)
    call mult_sq(g_,dx,LocalVa,8)
    do i=1,8
        g(i)=-g(i)*(sigma-sigma_n)-g_(i)*sigma
    enddo
  end
!*************************************************************************
subroutine local_dy(g,dy,h1,h2,LocalVn,LocalVa,sigma,sigma_n)
    dimension dy(8,8),LocalVn(8),LocalVa(8),g(8),g_(8)
    real*8 dy,LocalVn,LocalVa,g,g_,sigma,sigma_n,h1,h2

    dy(1,1)=-4.d0
    dy(1,2)=-2.d0
    dy(1,3)=4.d0
    dy(1,4)=2.d0
    dy(1,5)=-2.d0
    dy(1,6)=-1.d0
    dy(1,7)=2.d0
    dy(1,8)=1.d0

    dy(2,1)=-2.d0
    dy(2,2)=-4.d0
    dy(2,3)=2.d0
    dy(2,4)=4.d0
    dy(2,5)=-1.d0
    dy(2,6)=-2.d0
    dy(2,7)=1.d0
    dy(2,8)=2.d0

    dy(3,1)=-4.d0
    dy(3,2)=-2.d0
    dy(3,3)=4.d0
    dy(3,4)=2.d0
    dy(3,5)=-2.d0
    dy(3,6)=-1.d0
    dy(3,7)=2.d0
    dy(3,8)=1.d0

    dy(4,1)=-2.d0
    dy(4,2)=-4.d0
    dy(4,3)=2.d0
    dy(4,4)=4.d0
    dy(4,5)=-1.d0
    dy(4,6)=-2.d0
    dy(4,7)=1.d0
    dy(4,8)=2.d0

    dy(5,1)=-2.d0
    dy(5,2)=-1.d0
    dy(5,3)=2.d0
    dy(5,4)=1.d0
    dy(5,5)=-4.d0
    dy(5,6)=-2.d0
    dy(5,7)=4.d0
    dy(5,8)=2.d0

    dy(6,1)=-1.d0
    dy(6,2)=-2.d0
    dy(6,3)=1.d0
    dy(6,4)=2.d0
    dy(6,5)=-2.d0
    dy(6,6)=-4.d0
    dy(6,7)=2.d0
    dy(6,8)=4.d0

    dy(7,1)=-2.d0
    dy(7,2)=-1.d0
    dy(7,3)=2.d0
    dy(7,4)=1.d0
    dy(7,5)=-4.d0
    dy(7,6)=-2.d0
    dy(7,7)=4.d0
    dy(7,8)=2.d0

    dy(8,1)=-1.d0
    dy(8,2)=-2.d0
    dy(8,3)=1.d0
    dy(8,4)=2.d0
    dy(8,5)=-2.d0
    dy(8,6)=-4.d0
    dy(8,7)=2.d0
    dy(8,8)=4.d0

    do i=1,8
        do j=1,8
            dy(i,j)=dy(i,j)/72.d0*h1*h2
        enddo
    enddo

    call mult_sq(g,dy,LocalVn,8)
    call mult_sq(g_,dy,LocalVa,8)
    do i=1,8
        g(i)=-g(i)*(sigma-sigma_n)-g_(i)*sigma
    enddo
end
!*************************************************************************
subroutine local_dz(g,dz,h1,h2,LocalVn,LocalVa,sigma,sigma_n)
    dimension dz(8,8),LocalVn(8),LocalVa(8),g(8),g_(8)
    real*8 dz,LocalVn,LocalVa,g,g_,sigma,sigma_n,h1,h2

    dz(1,1)=-4.d0
    dz(1,2)=-2.d0
    dz(1,3)=-2.d0
    dz(1,4)=-1.d0
    dz(1,5)=4.d0
    dz(1,6)=2.d0
    dz(1,7)=2.d0
    dz(1,8)=1.d0

    dz(2,1)=-2.d0
    dz(2,2)=-4.d0
    dz(2,3)=-1.d0
    dz(2,4)=-2.d0
    dz(2,5)=2.d0
    dz(2,6)=4.d0
    dz(2,7)=1.d0
    dz(2,8)=2.d0

    dz(3,1)=-2.d0
    dz(3,2)=-1.d0
    dz(3,3)=-4.d0
    dz(3,4)=-2.d0
    dz(3,5)=2.d0
    dz(3,6)=1.d0
    dz(3,7)=4.d0
    dz(3,8)=2.d0

    dz(4,1)=-1.d0
    dz(4,2)=-2.d0
    dz(4,3)=-2.d0
    dz(4,4)=-4.d0
    dz(4,5)=1.d0
    dz(4,6)=2.d0
    dz(4,7)=2.d0
    dz(4,8)=4.d0

    dz(5,1)=-4.d0
    dz(5,2)=-2.d0
    dz(5,3)=-2.d0
    dz(5,4)=-1.d0
    dz(5,5)=4.d0
    dz(5,6)=2.d0
    dz(5,7)=2.d0
    dz(5,8)=1.d0

    dz(6,1)=-2.d0
    dz(6,2)=-4.d0
    dz(6,3)=-1.d0
    dz(6,4)=-2.d0
    dz(6,5)=2.d0
    dz(6,6)=4.d0
    dz(6,7)=1.d0
    dz(6,8)=2.d0

    dz(7,1)=-2.d0
    dz(7,2)=-1.d0
    dz(7,3)=-4.d0
    dz(7,4)=-2.d0
    dz(7,5)=2.d0
    dz(7,6)=1.d0
    dz(7,7)=4.d0
    dz(7,8)=2.d0

    dz(8,1)=-1.d0
    dz(8,2)=-2.d0
    dz(8,3)=-2.d0
    dz(8,4)=-4.d0
    dz(8,5)=1.d0
    dz(8,6)=2.d0
    dz(8,7)=2.d0
    dz(8,8)=4.d0

    do i=1,8
        do j=1,8
            dz(i,j)=dz(i,j)/72.d0 *h1*h2
    enddo
    enddo

    call mult_sq(g,dz,LocalVn,8)
    call mult_sq(g_,dz,LocalVa,8)
    do i=1,8
        g(i)=-g(i)*(sigma-sigma_n)-g_(i)*sigma
    enddo
end
!*************************************************************************
function mu_3d(i)
    mu_3d=mod(i-1,2)+1
end
!*************************************************************************
function nu_3d(i)
    nu_3d=mod((i-1)/2,2)+1
end
!*************************************************************************
function ku_3d(i)
    ku_3d=(i-1)/4+1
end
!*************************************************************************
!Adding local matrix & vector to global
subroutine AddLocal(di,ig,jg,ggl,ggu,pr,A,F,L,nEdge,mEdge,kEdge,igT,jgT,ggT,nodesTermCount,mTerm,nDim)
    dimension di(nEdge),ig(nEdge+1),jg(mEdge),ggl(mEdge),ggu(mEdge),A(kEdge,kEdge),&
    L(kEdge),F(kEdge,nDim),pr(nEdge,nDim),igT(nodesTermCount+1),jgT(mTerm),ggT(mTerm)
    real*8 di,ggl,ggu,A,F,pr,ggT,component

    do i=1,kEdge
        do j=1,kEdge
!           adding local matrix
            if(L(i).gt.nEdge.and.L(j).le.nEdge)then
                do k=igT(L(i)-nEdge),igT(L(i)-nEdge+1)-1
                    component=A(i,j)*ggT(k)
                    call addComponent(di,ig,jg,ggl,ggu,jgT(k),L(j),component,nEdge,mEdge)
                enddo
            else
                if(L(i).le.nEdge.and.L(j).gt.nEdge)then
                    do k=igT(L(j)-nEdge),igT(L(j)-nEdge+1)-1
                        component=A(i,j)*ggT(k)
                        call addComponent(di,ig,jg,ggl,ggu,L(i),jgT(k),component,nEdge,mEdge)
                    enddo
                else
                    if(L(i).gt.nEdge.and.L(j).gt.nEdge)then
                        do k=igT(L(i)-nEdge),igT(L(i)-nEdge+1)-1
                            do mu=igT(L(j)-nEdge),igT(L(j)-nEdge+1)-1
                                component=A(i,j)*ggT(k)*ggT(mu)
                                call addComponent(di,ig,jg,ggl,ggu,jgT(k),jgT(mu),component,nEdge,mEdge)
                            enddo
                        enddo
                    else
                        call addComponent(di,ig,jg,ggl,ggu,L(i),L(j),A(i,j),nEdge,mEdge)
                    endif
                endif
            endif
        enddo
!       adding local vector
        do in=1,nDim
            if(L(i).le.nEdge) then
                pr(L(i),in)=pr(L(i),in)+F(i,in)
            else
                do k=igT(L(i)-nEdge),igT(L(i)-nEdge+1)-1
                    pr(jgT(k),in)=pr(jgT(k),in)+F(i,in)*ggT(k)
                enddo
            endif
        enddo
    enddo
end
!*************************************************************************
subroutine addComponent(di,ig,jg,ggl,ggu,i,j,component,nEdge,mEdge)
    dimension di(nEdge),ig(nEdge+1),jg(mEdge),ggl(mEdge),ggu(mEdge)
    real*8 di,ggl,ggu,component
    if(i.eq.j) then
         di(i)=di(i)+component
    else
        if(i>j) then
            ibeg=ig(i)
            do while(jg(ibeg).ne.j)
                ibeg=ibeg+1
            enddo
            ggl(ibeg)=ggl(ibeg)+component
        else
            ibeg=ig(j)
            do while(jg(ibeg).ne.i)
                ibeg=ibeg+1
            enddo
            ggu(ibeg)=ggu(ibeg)+component
        endif
    endif
end
!*************************************************************************
subroutine local_bc_r(b,c,r1,r2)
dimension b(2,2),c(2,2)
real*8 b,c,r1,r2,h
    h=r2-r1

    b(1,1)=(r1+h/2.d0)/h
    b(1,2)=-(r1+h/2.d0)/h
    b(2,1)=-(r1+h/2.d0)/h
    b(2,2)=(r1+h/2.d0)/h

    c(1,1)=h/6.d0*(r1*2.d0+h/2.d0*1.d0)
    c(1,2)=h/6.d0*(r1*1.d0+h/2.d0*1.d0)
    c(2,1)=h/6.d0*(r1*1.d0+h/2.d0*1.d0)
    c(2,2)=h/6.d0*(r1*2.d0+h/2.d0*3.d0)
end
!*************************************************************************
subroutine local_bc_z(b,c,z1,z2)
    dimension b(2,2),c(2,2)
    real*8 b,c,z1,z2,h
    h=z2-z1
    c(1,1)=h/6.d0*2.d0
    c(1,2)=h/6.d0*1.d0
    c(2,1)=h/6.d0*1.d0
    c(2,2)=h/6.d0*2.d0

    b(1,1)=1.d0/h
    b(1,2)=-1.d0/h
    b(2,1)=-1.d0/h
    b(2,2)=1.d0/h
end
!*************************************************************************
function mu(i)
    mu=mod(i-1,2)+1
end
!*************************************************************************
function nu(i)
    nu=(i-1)/2+1
end
!*************************************************************************
subroutine v3d(vertex3D,nVer3D,vertex2D,nVer2D,vnRZ,vnXYZ)
    common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimEarth/nodesNotTermEarth,mEarth
    integer elemsCount2D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,&
    edgeTerm2d,searchElement2D,typeTask,calcWK,calcVP
    real*8 vertex2D,vertex3D,vnRZ,vnXYZ
    real*8 eps2d,eps3d,eps3dNS,radiusRED,x,y,z,r
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),vnRZ(nodesCount2D),&
    vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),vnXYZ(nodesCount3D)

    j=1
    do i=1,nodesNotTermEarth
        x=vertex3D(1,i)
        y=vertex3D(2,i)
        z=vertex3D(3,i)
        r=dsqrt(x**2+y**2)
        nEl=searchElement2D(vertex2D,nVer2D,r,z,nodesCount2D,elemsCount2D)
        if(nEl.ne.0)then
            call e2d(r,z,vertex2D,nVer2D,vnRZ,vnXYZ(j),1,nodesCount2D,elemsCount2D,nEl)
        else
            vnXYZ(j)=0.d0
        endif
        j=j+1
    enddo
end
!*************************************************************************
subroutine calculateA0(aX,aY,aZ,a,vertex3D,nodesforedges)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    real*8 vertex3D,a,aX,aY,aZ,axis,midValue,length
    integer elemsCount3D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,edgesBoundCount3D
    dimension vertex3D(3,nodesCount3D),nodesforedges(3,edgesCount3D),a(edgesCount3D),&
    aX(nodesCount3D),aY(nodesCount3D),aZ(nodesCount3D),axis(3),midValue(3)

    do i=1,edgesNotTermCount3D
        k1=nodesforedges(1,i)
        k2=nodesforedges(2,i)

        do j=1,3
            axis(j) = vertex3D(j,k2) - vertex3D(j,k1)
        enddo
        midValue(1) = (aX(k1) + aX(k2))/2.d0
        midValue(2) = (aY(k1) + aY(k2))/2.d0
        midValue(3) = (aZ(k1) + aZ(k2))/2.d0
        call scalar(axis,axis,3,length)
        length = dsqrt(length)
        call multVectorByNumberDouble(axis,1.d0/length,axis,3)
        call scalar(midValue,axis,3,a(i))
    enddo
end
!*************************************************************************
real*8 function ValueInPoint3D(v,vertex,nVer,x,y,z,num)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    real*8 vertex,v,x,y,z,x1,y1,z1,x2,y2,z2,functionX,functionY,functionZ,hx,hy,hz
    integer elemsCount3D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,edgesBoundCount3D
    dimension vertex(3,nodesCount3D),v(8),nVer(14,elemsCount3D)

    x1 = vertex(1,nVer(1,num))
    x2 = vertex(1,nVer(2,num))
    y1 = vertex(2,nVer(1,num))
    y2 = vertex(2,nVer(3,num))
    z1 = vertex(3,nVer(1,num))
    z2 = vertex(3,nVer(5,num))

    hx=x2-x1
    hy=y2-y1
    hz=z2-z1
    ValueInPoint3D=0.d0
    do i=1,8
        ValueInPoint3D=ValueInPoint3D+functionX(x1,x2,hx,x,i)*&
                                      functionY(y1,y2,hy,y,i)*&
                                      functionZ(z1,z2,hz,z,i)*&
                                      v(i)
    enddo
end
!*************************************************************************
!1-st bound condition
subroutine firstCondition(d,f,l1,n,nBound,nDim)
    dimension d(n),f(n,nDim),l1(nBound)
    real*8 d,f,p
    p=1e+30
    do j=1,nDim
        do i=1,nBound
            if(l1(i).le.n)then
                d(l1(i))=p
                f(l1(i),j)=0.d0
            endif
        enddo
    enddo
end
!*************************************************************************
!2-nd bound condition
    subroutine secondConditionForVn(f,i1,i2,i3,i4)
        common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
        m2D,mTerm2d,nodesNotTermCount2D
        dimension f(nodesNotTermCount2D)
        real*8 f
        !condition on Z (first source)
        f(i1)=-1.d0/(4.d0*3.1415926d0)
        f(i2)=-1.d0/(4.d0*3.1415926d0)
        !condition on Z (second source)
        f(i3)=1.d0/(4.d0*3.1415926d0)
        f(i4)=1.d0/(4.d0*3.1415926d0)
    end
    !*************************************************************************
!*************************************************************************
!*************************************************************************
